Ab-initio Phonon Calculations for the layered compound TiOCl 
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We present first-principles frozen-phonon calculations for the three Raman-active Ag modes in the spin- 1/2 
layered TiOCl system within two different well-known approaches: the local density approximation (LDA) 
and the so-called LDA+U approximation. We observe that the inclusion of electron correlation in a mean- 
field level as implemented in the LDA+U leads to a better overall agreement with experimental results. We 
also discuss the implications of the two approaches on the physics of TiOCl. 

PACS numbers: 75.30.Gw, 75.10.Jm, 78.30.-j 



Introduction - The layered quantum spin system TiOCl 
has been recently a subject of debate due to its anoma- 
lous properties at moderate to low temperatures. Suscep- 
tibility measurementsi show a kink at Tc2=94 K and an 
exponential drop at Tci=66 K indicating the opening of 
a spin gap which has been interpreted as a spin-Pcicrls 
phase transition!.. LDA+[/ calculations performed by 
these authors predicted the system to behave as a spin- 
chain along the b axis. X-ray diffraction measurements 
below Tci=66 K report^ the presence of superlattice re- 
flections at (h,k -\- 1/2, Z) which denotes a doubhng of 
the unit cell along the h axis and has been interpreted as 
a confirmation that the system undergoes a spin-Peierls 
phase transition below Td . Still unresolved is the nature 
of the kink at Tc2- The observation in Raman and in- 
frared spectroscopy of phonon anomalies^ above Tc2 as 
well as temperature-dependent g-factors and linewidths 
in ESR^ led to the proposal that the system is subject to 
competing lattice, spin, and orbital degrees of freedom. 
Recent heat capacity measurements^ seem to strongly 
support the idea of existence of possible lattice and/or 
orbital fluctuations above Tc2- 

In view of the above discussion, it is of great impor- 
tance to understand the behavior of the phonons in this 
system as well as the possibility of coupling to orbital and 
spin degrees of freedom. In a previous communication^ 
we presented first-principles results on the electronic 
properties of a few slightly distorted structures for TiOCl 
according to the various Ag-Raman active phonon modes 
expected for the Pmmn space group assigned to TiOCl. 
The idea in refp^ was to investigate the orbital occupancy 
as a function of lattice distortion, which could simulate 
the phonon-orbital interaction. The phonon modes and 
the corresponding eigenvectors in ref^ were obtained by 
considering a shell modeli. 

In the present work, we abandon the parametrized pro- 
cedure implicit in the shell model and we calculate the 
Raman-active Ag phonon modes in TiOCl within a first 
principles frozen phonon approach. From group theoret- 
ical analysis, the allowed Ag modes in the Pmmn sym- 
metry define ion displacements along the c-axis of the 
crystal. One important result of our calculations is that 
the consideration of correlations in this system is of cru- 
cial importance in order to get an accurate description of 
the phonons. 



First Principles Frozen Phonon approach.- In order to 
calculate phonon frequencies within an ab initio scheme, 
two different methods are generally employed: the energy 
surface method^ and the atomic force method^. Within 
both methods the harmonic potentials underlying the lat- 
tice vibrations are constructed by selecting frozcn-in dis- 
torted structures and then performing a Density Func- 
tional Theory (DFT) band-structure calculation. While 
in the first approach, energy values are needed to fit the 
energy surface to a bilinear form, in the second approach 
the atomic forces acting on the displaced ions are com- 
puted to fit the harmonic forces. 

In both methods the corresponding dynamical ma- 
trix is calculated according to the following procedure. 
Within the Born-Oppenheimer approximation, the elec- 
tronic energy hypersurface is represented as a function of 
the titanium, oxygen and chlorine coordinates along the 
c-axis. For small atomic displacements we approximate 
the surface by a Taylor expansion up to second order 
around its minimum (harmonic approximation). Since 
the minimum satisfies the equilibrium condition, the sur- 
face is described by the following bilinear form: 



= 1,2,3 
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where Zi are the c-axis coordinates of the ions (1 = 
titanium, 2 = oxygen, 3 = chlorine). The equation of 
motion for the i-th atom with mass Mi is then 
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Seeking for running- wave solutions of the form: Zi{t) 
_Jig-ii^t g^j,g jg£|. eigenvalue problem uj'^^i 
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where the matrix Dij — 



is the well- 



known dynamical matrix. 

Within the energy surface method, Eq. is exploited 
to obtain the matrix Z3y, while within the atomic force 
method, the left-hand side of Eq. [21is used. 

In this work we will make use of both methods within 
two different approximated forms for the exchange- 
correlation potential in the DFT scheme, namely, the 
Local Density Approximation (LDA)i& and the so-called 
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LDA+U approximationii. For the LDA+C/ calcula- 
tions, we considered the implementation, proposed 
by M.T.Czyzyk and G.A.Sawatzkjii^ (named in the liter- 
ature " AMF" , Around Mean Field) which differs with re- 
spect to the original Anisimov et al. energy functionaU^ 
in that it takes into account spin degrees of freedom ex- 
plicitly for all electrons in the system. Notably, the av- 
erage occupancy of one d orbital — 2(21+1) Smo- ^-mcr 
is now replaced by the spin dependent counterpart — 
2ITT '^"lo- and the starting energy functional is now 
the spin polarized LSDA functional, 

m.m' ,(7 

+\ E (c^-^)(«™.-«°)K.'.-n°). (3) 

m^m' .(7 

Calculations have been performed using the full- 
potential linearized augmented plane-wave code 
WIEN2k^^. While the atomic force method is more 
efficient to obtain the dynamical matrix, a successful 
implementation of this method within LDA-I-C/ is not yet 
reliable. Here we will present results on the structural 
properties of TiOCl within LDA considering the atomic 
force method and within LDA+J7 by means of the total 
energy method. 

Density Functional Calculations - In our calculations, 
the expansion of the wave functions included 1207 
LAPWe^s {RKmax = 7, including 12 Local Orbitals) 
and the muffin tin radii were chosen to be 1.7 a.u. for 
Titanium, 1.5 a.u. for Oxygen and 2.00 a.u. for Chlo- 
rine. Expansion in spherical harmonics for the radial 
wave functions were taken up to I = 10. Charge den- 
sities and potentials were represented by spherical har- 
monics up to L = 6 whereas in the interstitial region they 
were expanded in a Fourier series with 4500 stars. For 
Brillouin-zone (BZ) integrations a 500 k-point mesh was 
used yielding 60 k points in the irreducible wedge and 
use of the modified tetrahedron method was madei^. 

Concerning the LDA-I-C/ parameters U and J we have 
used the values of 4 eV and 1 eV respectively. Here we 
have considered a larger value of U than in the calcula- 
tion in refi&(U=3.3eV) in line with a Hartree-Fock calcu- 
lation of on-site Coulomb interaction for transition metal 
oxides-i. 

For a consistent phonon calculation we need to have a 
well defined ah initio equilibrium structure which corre- 
sponds to the minimum of the Taylor expansion Eq. ^ 
Therefore we analyze the structure of TiOCl by perform- 
ing an atomic position optimization and (within LDA) a 
unit cell volume optimization. In the volume optimiza- 
tion, the a : b : c ratios were kept constant and equal 
to the experimental values. For each total energy cal- 
culation we allow the ions to relax to their zero-force 
positions along the c-axis by simulating the dynamics of 
the damped Newton method. For the time evolution of 
the three atomic coordinates, the method is based on the 
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FIG. 1: Volume optimization for TiOCl within LDA with the 
ratio a : b : c fixed to the experimental one. 



TABLE I: Experimental and theoretical structure parameters 
(zTi, zo, zoj are the c-axis coordinates in units of the c-axis 
length and the distances Ti-0,Ti-Cl are given in A). 







zo 


ZCi 


Ti-0 


Ti-O-Ti 


Ti-Cl 


Exp. 


0.1194 


0.0551 


0.3318 


1.96 


150° 


2.40 


LDA 


0.0846 


0.0724 


0.2999 


1.89 


174° 


2.41 


LDA+U 


0.1133 


0.0523 


0.3223 


1.96 


151° 


2.38 



finite difference equation z^^"^ = zj + rj{zl — zl^^) + SF,l , 
where zf and F"^ are the coordinate and the force on the 
i-th atom at time step t. Damping and speed of motion 
are controlled by the two parameters rj and S respectively. 

In Fig. n the total energy is displayed versus the vol- 
ume variation with respect to the experimental reference 
value. The continuous line represents a fit to the Mur- 
naghan equation of statei^ from which we extract a bulk 
modulus value of 59.5 GPa. The curve shows that the 
optimal volume is in good agreement with the experi- 
mental one, in contrast with results obtained for high-T^ 
compounds^. 

In Table we show the equilibrium fractional c-axis 
coordinates, main distances and bonding angle of the 
three atoms obtained for the optimized lattice param- 
eters. The most important difference with respect to the 
experimental structure is the shape of the Ti-0 chain 
along the a-axis which in LDA has become almost a lin- 
ear chain (as shown also by the bonding angle Ti-O-Ti 
and the contraction of the Ti-0 distance). 

Recalling that the local Titanium coordinate frame has 
the z axis parallel to a and the x and y axes rotated by 
45° respect to b and c (see Fig. |5J, LDA provides a t2g 
ground state where the dxy and dxz , dyz contributions are 
comparable in weight (0.2 and 0.27 electrons/eV respec- 
tively). On the other hand, the ab initio calculations^ 
using the experimental structure indicate that the dxy 
weight is dominant . The approximate degeneracy of the 
three t2g orbitals is clearly shown in Fig.|21where the elec- 
tron density surface in the energy range -1 eV to Ep- for 
° 3 

an isovalue of O.le/A is drawn and refiects the typical 




FIG. 2: Electron density surface for an isovalue of O.leA for 
the LDA equilibrium structure (also shown the global and the 
local systems of reference). 



t2g spatial symmetry with the minima of the isodensity 
along the octahedron axis. The shortening of the Ti-0 
distance causes the local octahedron to be more regular 
and the shape of the surface around the titanium atom 
to be an equal admixture of dxy, dxz and dyz orbitals. 

Turning back to Tabled we observe that while the Ti- 
Cl distance is in good agreement with the experimental 
one, the Ti-0 distance is badly described by the LDA. We 
ascribe this failure of LDA to the improper description 
of electron correlations in titanium d-states and we show 
that inclusion of d-orbital correlations on a mean-field 
level (LDA+[/) causes the LDA equilibrium structure to 
be unstable with respect to the experimental one. 

In Fig. |31 we present the energy hypersurface in the 
fractional coordinate space (zTi,zo,zci) along the line 
connecting the LDA equilibrium position and the exper- 
imental one which is represented approximately by the 
vector (z,-z/2,z). 

It is clear that the inclusion of electron correlations 
distorts the LDA energy line in a sort of "double- well" 
curve whose minimum is now very close to the experimen- 
tal one (see Fig. 13 dotted line). Negative corrections to 
LDA total energy are expected for partially orbitally po- 
larized LDA groundstates, the amount of which depends 
on the correlation strength and on the spatial extension 
of the d orbitals (namely, the U value and the titanium 
muffin radius). In fact, as shown in Fig|21 the two energy 
curves match at the LDA equilibrium value of z where 
an orbitally non-polarized ground state is found. 

As a final step in the structure determination of TiOCl 
we want to find out the equilibrium structure according 
to the LDA+C/ functional. To determine the global min- 
imum of the energy hypersurface we adopt a trial and 
error approach by performing total energy calculations of 
all possible distorted structures around the experimental 
one and obtaining the minimum by a fitting procedure. 

As shown in Table the experimental distances and 
bonding angles are well reproduced by the LDA-|-[/ 
method; the disagreement in this case concerns the ab- 
solute fractional coordinates and it results, within ev- 
ery bilayer, in a small reduction of the distance between 
the two Ti-0 layers and accordingly between the outer 
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FIG. 3: LDA (solid line) and LDA+U (dashed line) energy 
curves along the line (z,-z/2,z) in the space (zTi,zo,zcO which 
connects the LDA equilibrium structure (taken as reference) 
to the experimental one. 
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FIG. 4: Electron density surface for an isovalue of O.le/A 
for the LDA+f/ equilibrium structure. 



CI ones. In Fig. 0] we show electron density plot of the 

" 3 

LDA-|-t/ equilibrium structure at the isosurface O.le/A . 

Preliminary calculations with the Generalized Gradi- 
ent Approximation^^ (GGA) method show that the op- 
timal equilibrium structure is improved with respect to 
the LDA results but it is still far from the experimental 
one. 

Phonons within LDA and LDA+U .- The phonon fre- 
quencies within the LDA approximation were calculated 
by considering the forces acting on the three atoms when 
they were displaced with respect to the equilibrium po- 
sitions in all possible fashions (in-phase and out-of-phase 
displacements of one, two and three atoms). The force 
values were fitted to linear polynomials in the c-axis co- 
ordinates and were determined within the uncertainty of 
1 mRy/a.u.. 

Table displays the LDA frequencies and the cor- 
responding relative atomic displacements in comparison 
with the Raman frequencies^. The lowest frequency 
mode is a titanium-chlorine in-phase and oxygen out- 
of-phase mode and it shows a considerable disagreement 
with respect to the experimental value. The second mode 
is a titanium-oxygen in-phase and chlorine out-of-phase 
mode and the third mode is mainly given by oxygen with 
small in-phase contribution of chlorine and almost neg- 
ligible out-of-phase contribution from titanium. To test 
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TABLE II: LDA and LDA+17 phonon frequencies compared 
to experimental data and mode assignment (relative atomic 
displacement, i.e. components of the eigenvectors divided by 
\JMi and normalized to the maximum amplitude). 







/3 


7 


c- 



Expt. 


froq. (cm ^ ) 


203 


365 


430 


LDA 


frcq.(cm^"'^) 


154 


351 


424 


(mode assign.) 


Ti 


0.909 


0.822 


-0.025 




O 


-0.157 


0.311 


1.000 




CI 


1.000 


-1.000 


0.104 


LDA+f/ 


freq.(cm^"'^) 


192 


336 


407 


(mode assign.) 


Ti 


0.557 


0.307 


-0.775 




O 


-0.048 


1.000 


1.000 




CI 


1.000 


-0.213 


0.614 



the accuracy of the fitting polynomials we crosschecked 
a posteriori the frequency values by displacing the atoms 
along the eigenvector directions finding the larger devia- 
tion (~ 10%) for the lowest mode. 

Within the LDA-t-J7 method we also calculated phonon 
frequencies and modes by means of a quadratic fitting of 
total energy surface. The frequencies and assignments 
are displayed in Table ^] The numerical error of the 
calculation is 10^'* Ry and affects mainly the second 
frequency value. The error due to the fitting function 
proved to be negligible confirming the harmonic character 
of all the three vibrations. In LDA+U the first frequency 
becomes considerably improved with respect to the LDA 
value and in good agreement to the experimental mode. 
The second mode upon inclusion of the numerical error 
is found in good agreement with the experiment and the 
third frequency disagrees by about 4%. 

Regarding the mode assignments in the two ap- 



proaches, both LDA and LDA-|-[/ calculations agree as 
far as the in-phase and out-of-phase features are con- 
cerned. From a quantitative point of view instead, 
LDA+C/ associates in contrast to LDA, a larger ampli- 
tude to chlorine and smaller to oxygen in the first mode; 
in the second mode the principal elongation is the oxygen 
one with smaller contributions from titanium and chlo- 
rine while the latter play the major role in LDA and fi- 
nally the third mode is characterized by roughly balanced 
amplitudes for all the three atoms in LDA-I-C/ while in 
LDA the mode is almost pure oxygen oscillation. Calcu- 
lation of the Raman cross sections for the three modes 
would bring to evidence the main differences between the 
LDA and LDA-|-?7 atomic displacements which can be 
tested against the experimental Raman intensities. 

Conclusions.- Summarizing, we have presented DFT 
calculations for the Ag Raman phonon modes in TiOCl 
within LDA and LDA+C/. Within LDA we obtain that 
the optimized equilibrium structure deviates consider- 
ably respect to the experimental one while considera- 
tion of electron correlation within the hDA+U approach 
brings the DFT minimum very near to the experimental 
data (importance of correlation effects has been already 
pointed out by P. Labeguerie et alM- and by M. Mer- 
awa et al.^). In both approaches, by calculation of the 
dynamical matrix, we obtain the three Ag modes in the 
Pmmn symmetry. While both LDA and LDA-|-[/ show 
an overall qualitative agreement with the results from 
Raman scattering experiments, quantitatively LDA-I-U 
produces for the lower mode the best agreement with 
the experiment while LDA performs better for the sec- 
ond and third mode. 
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